2. Meshing from Neper geometries

Neper is a free / open source software package for polycrystal generation and meshing, see https://neper.info

To run this notebook you first have to install neper.

from netgen.occ import *
from ngsolve import *
from ngsolve.webgui import Draw
import subprocess, random

n = 100
subprocess.run(["neper", "-T", "-n", str(n), "-format", "tess"])

========================    N   e   p   e   r    =======================
Info   : A software package for polycrystal generation and meshing.
Info   : Version 4.9.1-6
Info   : Built with: gsl|muparser|opengjk|openmp|nlopt|libscotch (full)
Info   : Running on 8 threads.
Info   : <https://neper.info>
Info   : Copyright (C) 2003-2022, and GNU GPL'd, by Romain Quey.
Info   : No initialization file found (`/Users/joachim/.neperrc').
Info   : ---------------------------------------------------------------
Info   : MODULE  -T loaded with arguments:
Info   : [ini file] (none)
Info   : [com line] -n 100 -format tess
Info   : ---------------------------------------------------------------
Info   : Reading input data...
Info   : Creating domain...
Info   : Creating tessellation...
Info   :   - Setting seeds... 
Info   :   - Running tessellation...
Info   : Generating crystal orientations...
Info   : Writing results...
Info   :     [o] Writing file `n100-id1.tess'...
Info   :     [o] Wrote file `n100-id1.tess'.
Info   : Elapsed time: 0.017 secs.
========================================================================

CompletedProcess(args=['neper', '-T', '-n', '100', '-format', 'tess'], returncode=0)

verts = {}
edges = {}
faces = {} 
solids = []

f = open("n" + str(n) + "-id1.tess")

def randomcol():
    r = random.uniform(0, 1)
    g = random.uniform(0, 1)
    b = random.uniform(0, 1)
    return (r,g,b)


while True:
    line = f.readline()
    if not line: break

    if line.split()[0] == "**vertex":
        num = int(f.readline())
        print ("found", num, "vertices")

        for i in range(num):
            line = f.readline()            
            nr,x,y,z,hhh = line.split()
            verts[int(nr)] = Vertex(Pnt(float(x), float(y), float(z)))

    if line.split()[0] == "**edge":
        num = int(f.readline())
        print ("found", num, "edges")

        for i in range(num):
            line = f.readline()            
            nr,i1,i2,hhh = line.split()
            edge = Edge(verts[int(i1)], verts[int(i2)])
            edges[int(nr)] = edge
            edges[-int(nr)] = edge.Reversed()
            
    if line.split()[0] == "**face":
        num = int(f.readline())
        print ("found", num, "faces")

        for i in range(num):
            l1 = f.readline()            
            l2 = f.readline()            
            l3 = f.readline()            
            l4 = f.readline()
            nr,*_ = l1.split()
            nume,*fedges = l2.split()
            
            face = Face(Wire( [edges[int(enr)]  for enr in fedges] ))
            faces[int(nr)] = face
            faces[-int(nr)] = face.Reversed()

            
    if line.split()[0] == "**polyhedron":
        num = int(f.readline())
        print ("found", num, "polyhedra")

        for i in range(num):
            nr,nrf,*polyfaces = f.readline().split()
            solids.append (Solid(Glue( [faces[int(fnr)] for fnr in polyfaces ] )))
            solids[-1].faces.col = randomcol()                        
            
shape = Glue(solids)
geo = OCCGeometry(shape)

Draw (shape);

found 575 vertices
found 1146 edges
found 672 faces
found 100 polyhedra

with TaskManager(pajetrace=10**9):
    mesh = geo.GenerateMesh(maxh=0.1)

print ("ne =", mesh.ne)

ne = 60279

Draw (mesh);